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ABSTRACT 

A new highly accurate algorithm for the solution of the Falkner-Skan equation of 
boundary layer theory is presented. The algorithm, based on a Maclaurin series 
representation, finds its coefficients from recurrence. In addition, Wynn-epsilon 
convergence acceleration and continuous analytical continuation enable an 
accurate evaluation. The most accurate skin friction coefficients (shooting angle) 
to date are presented along with comparisons to past and present values found in 
the literature. The algorithm, coded in FORTRAN, uses neither enhanced 
precision arithmetic beyond quadruple precision nor computer algebra to achieve 
results in a timely fashion. 

Key Words: Falkner-Skan flow; Blasius flow; Wynn-epsilon acceleration; 
Romberg acceleration; Continuous analytical continuation 

1. INTRODUCTION 

Every since its first appearance in the literature in 1908 [1], the Blasius equation 
describing viscous flow over a flat plate has fascinated physicists, engineers, 
mathematicians and numerical analysts alike. This ODE is rich in physical, 
mathematical and numerical challenges. Two-dimensional flow over a fixed 
impenetrable surface creates a boundary layer as particles move more slowly near 
the surface than near the free stream. Because of its application to fluid flow, 
physicists and engineers have a keen interest in solving the Blasius equation and 
the related, but more general, Falkner-Skan (F-S) equation [2]. Since one can 
elegantly reduce these equations to one-dimensional non-linear ODEs through 
similarity arguments, mathematicians have found their fulfillment in uncovering 
the underlying symmetries and proving existence and (non-) uniqueness of its 
solutions. Unfortunately, a general analytical solution has not been forthcoming; 
however, for special cases of the F-S equation, several analytical solutions do exist. 
These prove most beneficial in verifying numerical algorithms. Numerical 
analysts, or as they are called by J.P Boyd, "arithmurgists" [3], have had a "field- 
day with these equations. They offer the mystery of nonlinearity — yet the 
simplicity of one-dimension — and the challenge of solving a boundary value 



problem through the determinism of an initial value problem. A host of numerical 
methods have emerged to solve these equations including, but not limited to, finite 
differences and finite elements, Adomian's polynomials, perturbation methods, 
homotopy, differential transformation, generalized Laguerre and Chebyshev 
polynomial expansions, variational iteration, the quasi-linear approximation and 
the diagonal Pade approximant. A list of nearly 150 references (some repeated) for 
these and other algorithms is found in Refs. 4 and 5. We also note that the 
majority of the numerical algorithms are based on Runge-Kutta ODE solvers, and 
therefore have a definite "black box" quality. In spite of the enormous numerical 
effort however, a truly simple, yet numerically accurate and robust algorithm is 
still missing. Many, if not all, algorithms to this point seem rather delicate in that 
their iterative strategies must be carefully tuned to avoid numerical instability. For 
example, most schemes require the initial guess of the shooting angle to be 
relatively close to the converged result, which does not make for a robust 
algorithm. Judging from recent literature, the general lack of numerical agreement 
to a consistent five or more digits is indicative of the need for a reliable algorithm— 
the development of which we now address. 

The goal of the present investigation is to provide a robust algorithm for the 
solution of the F-S equation for physically relevant flows and for some additional 
flows found in the literature. We report highly accurate solutions to at least 10 or 
more digits without resorting to enhanced precision arithmetic or computer 
algebra. All calculations are programmed in FORTRAN double precision (DP) 
arithmetic with the exception of several quadruple precision (QP) demonstrations 
of claimed accuracy. 

We consider the F-S equation, 

rW^o/W/"W + ^(l-/'W 2 ) = 0, 7g[0,°o), (la) 

Physically, the F-S equation describes 2D flow over stationary impenetrable wedge 
surfaces of included angle fin, which limits to a flat plate, and the Blasius solution, 
as /? approaches zero. Eq(la) is subject to the following boundary conditions: 

/(°) = ° 

/'(0) = (lb) 
lim/'(i7) = l, 

7-»oo 



where the velocity profile goes asymptotically to unity. While 
/3^ a □ 37735 < P < 2 represents physical flow, with two solutions in the 

negative range of /3 for forward and reverse flow, solutions for larger values of /?, 
as well, found in the literature, are also considered. 

Curiously, as will be shown, we find that the most numerically accurate form of 
solution was originally proposed by Blasius himself in 1908 [1] as the Maclaurin 
series 

/M = I^7fV- (2) 

Moreover, all the elements of a particularly robust solution, exhibiting superior 
accuracy already exist in the literature. The relations between the elements along 
with some numerical common sense then provide the desired robust numerical 
algorithm we seek. 

In the following section, we derive the fundamental algorithm based on recurrence 
for the derivatives in the solution representation of Eq(2). Note that the recurrence 
was derived in Ref. 6 in an equivalent vector form, but the chosen implementation 
limited the method's full potential. Even with its coefficients known, direct 
evaluation of Eq(2) is problematic because of the finite radius of convergence of 
the series representation [7]. For this reason, as purposed in Ref. 8, we apply a 
Pade approximant. Unlike that suggested however, we seamlessly implement the 
diagonal approximant through Wynn-epsilon convergence acceleration. 
Unfortunately, convergence acceleration alone is insufficient to give the desired 
high accuracy. To achieve high accuracy therefore, continuous analytical 
continuation (CAC), applied in the manner suggested also in Ref. 6, is the remedy. 
Then, with confidence in our numerical evaluation, we are in position to 
implement the shooting method taking advantage of the physical features of the 
velocity profile as identified in Ref 9. We conclude with comparisons to literature 
values and, when convenient, add digits to establish new benchmark standards. 

Therefore as indicated, an accurate algorithm for resolution of the F-S equation has 
existed in the literature since the application of Wynn-epsilon convergence 
acceleration [10]. Its emergence however, required reliable computing, courage to 
embrace experimental mathematics and an extreme numerical desire. 



2. METHOD OF SOLUTION 

In this section, we derive the recurrence for the derivatives of Eq(2) and formulate 

the shooting method. 

2.1 Maclaurin series representation 

A recurrence relation for the derivatives in the Maclaurin series is found by 
applying the differential operator [6] 

d k 



dri k 

toEq(la). In the process, note that 

and 

to give 



(3) 



If 



f k) (n) = k\a t , 

then, with some algebra 



a k M = k (k-\) ^ ai (V)^ {*?)- A ( k ~ l ) a o frK-i (^)] "f + 

(4a) 

and the boundary conditions become 

a (o) = o 

a 1 (0) = 
lim a x (77) = 1. 



Thus, 



/w=z%(oy (5) 

is the desired evaluation. 

There are three issues to emphasize here. First, the dependent variables remain in 
their original domain and are not in a finite interval as found in many of the 
previous investigations. Such a transformation gives a false sense of security of 
improved accuracy because one has avoided the infinity of r\. With a 
transformation however comes the need to include an additional boundary 
condition, which unnecessarily complicates the shooting strategy by requiring an 
optimum search in a 2D parameter space (for the shooting angle and the additional 
condition). The second point is that by formulating the solution in terms of 
derivatives only, the ODE's nonlinearity is accommodated naturally through the 
(nonlinear) recurrence itself. There is no need to explicitly consider the 
nonlinearity of Eq(la) any further. Finally, we are not reliant on an ODE solver 
with its inherent truncation and roundoff errors. While marching is still required, it 
is more accurately accomplished by evaluating power series solutions over 
contiguous intervals. 

2.2 THE SHOOTING METHOD 

The shooting method is the preferred way to treat the F-S boundary value problem. 
Through a specification of an additional initial condition to replace the condition at 
infinity, the boundary value problem transforms into an equivalent iterative initial 



value problem. The initial condition is f"(0,a) = a , where a is the skin friction 
coefficient. To be equivalent, the shooting angle a must be determined such that 
lim/'(77;a) = l. Except where noted for a particular range of /?, the solution is 

t]— >oo 

assumed unique [12]. Now, the recurrence begins with 

a Q (0,a) = 

fl 1 (0,a) = (4b) 
a 2 (0,a) = a 12. 

In this way, Eq(5), with a assumed, gives / and its first two derivatives as 

00 

f(f];a) = ^a k (0,a}] k 

00 

r(/ 7 ;«)=iH(o,«y- 1 . (?) 

00 

f"(r ] ;a) = Zk(k-l)a k (0,ay J k -\ 

k=0 

where a prime refers to differentiation with respect to rj. Note the explicit 
dependence of a since a is to be varied. 

3. NUMERICAL IMPLEMENTATION OF THE SERIES SOLUTION 
REPRESENTATION 

We now consider how to best evaluate Eq(5) to a desired accuracy, which will be 
through convergence acceleration. 

3.1 Evaluation of the Maclaurin series by convergence acceleration 

It is well known that the radius of convergence for the series in Eq(5) is finite, thus 
limiting its usefulness. For this reason, we apply Wynn-epsilon convergence to 
accelerate the partial sums to convergence. 

First, consider the evaluation of a general infinite series 



00 

k=0 



which is, in the limit 



5(/7) = lim^(/7), 
if it exists, where the partial sums are 

«,(7) = iw'. 

k=0 

In this way, a sequence of approximations to the limit 

S = {5 / (; 7 ),/ = l,2,....}, 
is created and assumed to uniformly converge. 

To accelerate convergence, one applies the Wynn-epsilon (We) [11] convergence 
accelerator, which frequently results in a more rapidly converging sequence than 
the original. In this regard, we enter into the realm of experimental numerical 
methods since there is no proof that the F-S series solution lends itself to such an 
acceleration. There is some indication of success however, in that a diagonal Pade 
approximant gives a solution [8] and the We algorithm is just that [13]. 

The We accelerator takes the following form for the sequence S l (77); / = 0,1,.. .,L : 
£ { Ss i (r ] ),l = Q,...,L 

, k = 0,...,L ; l = Q,...,L-k-\. 
The recurrence forms a tableau 



o( /+1 )_,(0 
b k b k 



'0 

,(!) 

'0 



,(o) 
'1 

,(1) 



,(o) 

'2 
,(1) 



,(^-2) 



.(0) 
'L-l 

,« 
'L-l 



,(0) 



(8) 



where each element of an even column estimates the limit. Convergence comes 
from interrogation of the last term of even columns s\ L ~'\ i = Q,2,...,2[L/2] 



e We ~ 



£ {L-i) _ 8 {L-i-2) 



XL-i) 



<s, i = 2,...,2[L/2]. 



(9) 



Therefore, with Si appropriately defined, the limits of the sequences 



fi(?7,a) = Y J a k (0,aW 



k=0 



fi'(?7,a) = j:ka k (0,ay 7 k - 1 



k=0 



f l "(j 1 ,a)=? j k{k-\)a k (0,a) n 



k-2 



k=0 



(10) 



give the desired solution, provided all converge, which also has not been shown, 
but is obviously assumed. 

As a verification of convergence acceleration for f5 > , we first consider the 
Blasius solution [ (5 = ] since there is a known highly accurate value of a for this 
case. In particular, the shooting angle from an extensive Runge-Kutta calculation 
is 



a = 0.33205733621519630. 



(11) 



to 17 digits [7]. Then, assuming this value initiate the recurrence relation of 
Eq(4a) with Eqs(4b) gives Fig. 1 for the convergence of the velocity profile partial 




Fig. 1 . Sequence convergence for the Blasius solution for the velocity profile 
from Eq(5) without (a) and with (b) We acceleration. 




5 10 15 20 25 30 



1 

Fig. 2a. Velocity profile f'{rj,a) from Eq(5). 

sums fl without (a) and with (b) We acceleration. The effectiveness of We 

acceleration is quite apparent. As 77 approaches the radius of convergence, the 
original partial sums begins to oscillate while the We converged sequence remains 
steady and accurate. If we allow r\ to become larger than radius of convergence, 
Fig. 2 results for The original sequence is now entirely out of 

numerical control. However, the We accelerated sequence converges to 5 or 6 
places of the correct asymptotic value even when the original sums are take values 
of 10 14 or 10 15 . This is a remarkable demonstration, in general, of the power of 
convergence acceleration. 

However, eventually (for r\ > 8), even the We sequence fails to converge. 
Therefore, to overcome this difficulty, we propose continuous analytical 
continuation. 

3.2 Continuous analytical continuation (CAC) 

CAC [6] is the continuation of a Taylor series into the complex plane, or in this 
case along the real axis. This is accomplished by re-starting the series at successive 




-I -9. 
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Fig, 2b. Velocity profile f'(rj,a) from CAC. 
intervals, A?]. Thus, Eq(4a) remains valid but now the initial three derivatives are 

a l {rio' a ) = f'{ r l^ a ) (4b) 

which replace Eq(4b). % is the initial point of the interval, and 
f(r] ,a),f'(r] ,a),f"(r] ,a) are presumably available from the end of the 
previous interval. The desired quantities are therefore 

f(?] + Aj],a) = ^a k (j] ,a)(Aj]) 

k=0 



00 

f'(j] +Aj],a) = ^ka k (j] ,a)(A?]) 



(12) 



cc 



k-2 



f(j 1o +Aj 1 ,a) = ^k(k-\)a k (j] ,a)(Aj]) 



With CAC applied to the case of Fig 2a, Fig 2b results. Now the unaccelerated 
(original) partial sums converge up to tj = -45, but eventually fail. The We 
converged sequence however, maintains the correct value to 77 = 125 and beyond — 
again a remarkable demonstration of We acceleration. It would certainly be a 
challenge for a Runge-Kutta based method, or any other existing method for that 
matter, to reproduce this degree of accuracy. In particular, the challenge is for a 
numerical method, with no more than 16-digit precision and a mesh spacing of 1, 
to give the Blasius velocity profile at unity to eight digits of accuracy for 77 greater 
than 100. The accuracy of the evaluation of Eq(5) at the initial point of the interval 
via converged acceleration, is the key to the success of the evaluation. 

4. THE SHOOTING ALGORITHM 

We now formulate the conventional shooting algorithm. The concept is to assume 
the shooting angle a and iterate until the condition 



is satisfied to an acceptable accuracy. One of the difficulties of this approach is the 
sensitivity of the velocity profile to a. We develop the algorithm for P > and 
P < separately. 

4.1 For p > 

Consider Homann flow near the stagnation point of the boundary layer on a surface 
of revolution. For this case, P = 2, P — 1 , where a 1 0-place value of a is known 
to be [14] 




(13) 



a = 1.3119376939. 



(14) 



Figure 3a shows the velocity profile as a varies from 2 to 0. As we pass through 
the correct a [Eq(14)], the profile takes on a value of unity at large 77, which is 
more evident in the tight view in Fig 3b. This behavior seems typical of the general 




10 
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Fig 3. Variation of velocity profile with a for Homann flow: 
(a) Wide view (b) Tight view. 



case for J3 > and is true for all the physically relevant (and physical) J3 
investigated. The only differences occur for 77 usually under 10. In particular, we 
pass through the asymptotic velocity profile by decreasing a from above as has 
been shown numerically for positive /3 [9]. In addition, the asymptotic profile is 
approached from below without overshoot [9]. While this has not been shown for 
all positive values of /? , from this author's experience it seems true and therefore 
is assumed until circumstances dictate otherwise. The shooting algorithm therefore 
consists of the following procedure: 

a. Starting from a relatively large value of a as the initial guess, 
f'{rj) is evaluated by increasing 77 through steps of h from 
zero to rj m . 

b. If at some 77, /'(//) >1, then a is decreased and f'yl) evaluated 
until y (77) < 1 for some a. At this point, the asymptotic profile 
is bracketed. 

c. A new a is then determined by bisection and interrogation of the 
profile continues with subsequent a's found by bisection as f'(r/) 
oscillates about unity. 

d. If f'(ri) does not cross unity from below as r\ increases from 

zero to 77 m , then f"(rj) is checked for negativity. If negative, 

a is below its correct value and bisection again determines the 
next a. Experience indicates that r/ m = 20 is sufficient for /?> 0. 

e. Finally, when the estimate for a is approximately within an order 
of magnitude of the desired error, we apply the secant method 
rather than bisection. 

Table 1 gives results for the above procedure applied to Homann flow. The entry 

1 2 

for error of 10" confirms the value quoted by Zhang and Chen [14] over that of 
Asaithambi [6]. Also indicated is the number of (bisection) iterations for 
convergence. It should be noted that bisection is the preferred method over 
Newton-Raphson to find the correct a because of its reliability. Hence, the 
security of bisection is more highly valued than the faster convergence of Newton- 
Raphson. In this way, potential instability is avoided. In any case, to obtain the 
accuracy of the recent literature, less than 30 iterations are required. By all 
accounts, about 50 iterations are typical for 14-digit accuracy. 

As a further indication of the power of the convergence acceleration, a case, at an 



Table 1 



Variation of a for Homann Flow with the specified error 
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#Itr 
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10" 9 


1 31193769 


28 


io- 10 


1.311937690 


33 


10" 11 


1.3119376938 


34 


IO" 12 


1.31193769392 


40 


IO" 13 


1.31193769388 


39 


IO" 14 


1.3119376938798 


44 


io- 32 


1.3119376938798051354816461707 


104(QP) 



extreme accuracy of 10" (in Quad Precision arithmetic) is included. 29 digits are 
expected to be accurate (last entry) as verified by the Pohlhausen flow example 
considered next. 

Pohlhausen flow, for which j3 Q = O, (3 = 1 , is a rare instance of an analytical 

solution to the F-S profiles and for a, which is 21 S- As shown in Table 2, the 
proposed algorithm gives 1 5 correct places for a requested error of 10" 16 and 29 

Table 2 

Pohlhausen flow /3 = 0, (3 = 1 

e a 

IO" 16 1.154700538379252 
1 0" 32 1 . 1 5470053 837925 1 5290 1 829756 1 1 
Exact 1 . 1 5470053 83 7925 1 5290 1 829756 1 00 



(QP) for an error of 10" . Arguably, this example provides a most convincing 
demonstration of convergence acceleration coupled to continuous analytical 
continuation. 

Before addressing the standard cases found in the literature, we consider the 
Blasius solution in some detail. Table 3a gives highly accurate values of a and 
Table 3b gives the Blasius function and its derivatives to 10-digits. These values 
are in exact agreement with the 8-place results of in Ref. 15, which seems to be the 
most accurate published to date. Please note that the values in this reference were 
obtained from an integral formulation with Eq(5) at its core in 1981. In addition, 



Table 3a 

Blasius Flow /3 = 1, J3 = 



a #Itr 



Exact 0.33205733621519630 

s= 10" 16 0.332057336215195 52 
£=10" 32 0.33205733621519629893718006201 104 (QP) 



the asymptotic value 

lim [77 -/(;;)] 



is 1.7207876575205 in DP arithmetic agreeing to 14-digits with the value of Ref. 
3. In QP arithmetic, we find 1.720787657520502812, which agrees to all digits 
when rounded and adds several. The degree of agreement for a calculation without 
multiple precision arithmetic beyond QP is rather notable and should be the 
standard to which algorithms are held. 

All calculations for the above 3 figures and 3 tables required under 3 minutes on a 
Gateway 1.4 GHz laptop. 

4.2 For 

Figure 4 gives the variation of the velocity profile as a varies from 1 to for 
P = -0.1. As for positive /?, when a passes through its correct value [9], the 

velocity takes on an its asymptotic profile. However, unlike for positive /?, the 
velocity profile initially overshoots its asymptotic limit. Even with this difference, 
the shooting strategy outlined above remains valid since the approach is still from 
above. 



5. COMPARISONS OF RESULTS 

In this section, we compare the results of the present algorithm for the F-S solution 
with some of the most accurate results found in the literature. The physical flows, 
forward and reverse, are considered separately. 

5.1 Forward flows 

Table 4a lists converged a for positive /3 to 40. When rounded and compared to 
the corresponding results of Refs. 14, 16 and 17, we find all values in agreement to 
6-places except for the last digit as indicated in Table 4b. In addition, the free 
boundary rj^ layer estimate, defined here when the velocity profile becomes unity 



Table 3b 

Blasius Flow Profiles 



n 



f 



f" 



0.0E+00 

2.0E-01 

4.0E-01 

6.0E-01 

8.0E-01 

1.0E+00 

1.2E+00 

1.4E+00 

1.6E+00 

1.8E+00 

2.0E+00 

2.2E+00 

2.4E+00 

2.6E+00 

2.8E+00 

3.0E+00 

3.2E+00 

3.4E+00 

3.6E+00 

3.8E+00 

4.0E+00 

4.2E+00 

4.4E+00 

4.6E+00 

4.8E+00 

5.0E+00 

5.2E+00 

5.4E+00 

5.6E+00 

5.8E+00 

6.0E+00 

6.2E+00 

6.4E+00 

6.6E+00 

6.8E+00 

7.0E+00 

7.2E+00 

7.4E+00 

7.6E+00 

7.8E+00 

8.0E+00 

8.2E+00 

8.4E+00 

8.6E+00 

8.8E+00 



0.000000000E+00 
6.640999715E-03 
2.655988402E-02 
5.973463750E-02 
1.061082208E-01 
1.655717258E-01 
2.379487173E-01 
3.229815738E-01 
4.203207655E-01 
5.295180377E-01 
6.500243699E-01 
7.811933370E-01 
9.222901256E-01 
1.072505977E+00 
1.230977302E+00 
1.39680823 1E+00 
1.569094960E+00 
1.746950094E+00 
1. 929525 170E+00 
2.11 60298 17E+00 
2.30574641 8E+00 
2.498039663E+00 
2.692360938E+00 
2.888247990E+00 
3.085320655E+00 
3.283273665E+00 
3.481867612E+00 
3.680919063E+00 
3.880290678E+00 
4.079881939E+00 
4.279620923E+00 
4.479457297E+00 
4.679356615E+00 
4.87929581 1E+00 
5.079259772E+00 
5.27923881 1E+00 
5.479226847E+00 
5.679220147E+00 
5.8792 16466E+00 
6.07921448 1E+00 
6.27921343 1E+00 
6.479212887E+00 
6.679212609E+00 
6.879212471E+00 
7.079212403E+00 



0.000000000E+00 
6.6407792 10E-02 
1. 32764 1608E-01 
1.989372524E-01 
2.647091387E-01 
3.297800312E-01 
3.937761044E-01 
4.5626 17647E-01 
5.167567844E-01 
5.747581439E-01 
6.297657365E-01 
6.813 103772E-01 
7.289819351E-01 
7.72455021 1E-01 
8.115096232E-01 
8. 46044443 7E-01 
8.760814552E-01 
9.017612214E-01 
9.233296659E-01 
9.411179967E-01 
9.555182298E-01 
9.669570738E-01 
9.758708321E-01 
9.826835008E-01 
9.877895262E-01 
9.915419002E-01 
9.942455354E-01 
9.961553040E-01 
9.974777682E-01 
9.983754937E-01 
9.989728724E-01 
9.993625417E-01 
9.99611 70 17E-01 
9.997678702E-01 
9.998638190E-01 
9.999216041E-01 
9.999557173E-01 
9.999754577E-01 
9.999866551E-01 
9.999928812E-01 
9.999962745E-01 
9.999980875E-01 
9.999990369E-01 
9.999995242E-01 
9.999997695E-01 



3.320573362E-01 
3.319838371E-01 
3.314698442E-01 
3.300791276E-01 
3.273892701E-01 
3.230071 167E-01 
3.16589191 1E-01 
3.078653918E-01 
2.9666346 15E-01 
2.829310173E-01 
2.667515457E-01 
2.483509132E-01 
2.2809 17607E-01 
2.064546268E-01 
1.840065939E-01 
1.613603 195E-01 
1.391280556E-01 
1.178762461E-01 
9.808627878E-02 
8.012591814E-02 
6.423412109E-02 
5.051974749E-02 
3.897261085E-02 
2.948377201E-02 
2.1871 18635E-02 
1.590679869E-02 
1.134178897E-02 
7.9276598 15E-03 
5.431957680E-03 
3.648413667E-03 
2.402039844E-03 
1.550170691E-03 
9.806151 170E-04 
6.080442648E-04 
3.695625701E-04 
2.201689553E-04 
1.285698072E-04 
7.359298339E-05 
4.129031 111E-05 
2.270775 140E-05 
1.224092624E-05 
6.46797861 1E-06 
3.349939753E-06 
1.700667989E-06 
8.462841214E-07 



A) = 1,/?=-0.1 



2 4 6 8 10 12 14 

Fig.4. Variation of the velocity profile for /?= -0.1. 



Table 4a 

For < P 



p 


a 




#Itr 


40 


7.31478497433 


1.57-1.58 


49 


30 


6.33820862834 


1.79-1.80 


57 


20 


5.18071802491 


2.14-2.15 


50 


15 


4.49148689764 


2.42-2.43 


48 


10 


3.67523410111 


2.83-2.84 


53 


2 


1.68721816921 


4.46-4.47 


49 


1 


1.23258765682 


4.98-4.99 


53 


0.5 


0.927680039837 


5.37-5.38 


50 





0.469599988361 


6.07-6.08 


52 



Table 4b 

Last digits not in agreement 
inRefs. 14,16 and 17 
for /?s indicated 

14 1,15,30 

16 1 
,11 -0.15 

^<0{16 -0.18 

17 -0.15,-0.18 



^>0 



Table 5 

Forward flow for -0.198837735 < j3 < 






a 


TJoo 


Itr 


-0.10 


3.19269759843 


6.36-6.37 


48 


-0.15 


2.16361405647 


6.60-6.61 


48 


-0.18 


1.28636220596 


6.85-6.86 


49 


-0.1988 


5.218187884E-03 


7.32-7.33 


42 


-0.198837 


7.24675233 E-04 


7.32-7.33 


59 


-0.1988377 


1.58136616 E-04 


7.35-7.36 


91QP 


-0.198837735 


5.77016 E-06 


7.35-7.36 


95QP 



to 5 places, is also given in Table 4a. Finally, as a measure of the computational 
effort, the number of bisections for DP arithmetic is also provided. 

Table 5 gives results for fi^ < (5 < . For the first three entries, all values agree 
to the fifth place with Refs. 14, 16 and 17 except as again indicated in Table 4b. 
For the fourth entry, only the first 2 digits agree however. This value of ft is near 
the boundary between existence and non-existence (/Lin) of the solution and is 
particularly difficult numerically troubling. The value quoted in the table is 
expected to be accurate to all but the last digit. The next three entries are an 
attempt to more accurately define this boundary. For these /?'s, the determination 
of a requires QP arithmetic and further independent verification. 

The final comparison for forward flows is with a relatively obscure but very 
important data set. Apparently, M. Katagiri [18], in a symposium on space flight 
held at the University of Tokyo, published what must be considered the most 



accurate determination of a before the present investigation. Using a quasi-linear 
approximation in an integral formulation, Katagiri was able to compute values of a 
for -0.19</?<1 to 10 or 11 digits. Unfortunately, this work has been lost in 

antiquity, but its accuracy has been resurrected here in Table 6 for a with two 
additional digits as found from the present algorithm. All 40 values are in roundoff 
agreement with Katagiri 's results. The table required 2 min of computational time 
on a Gateway 1 .4 GHz laptop. This agreement makes the work of Katagiri that 
much more impressive as Katagiri published his results in 1986— during the era of 
primitive computing. In addition, had this reference been more widely known in 
the West, many of the recent numerical approaches would never have been 
attempted, as they are far inferior. It is only fitting to acknowledge Katagiri as the 
first numerical analyst to achieve truly high precision results for the solution of the 
F-S equation. 

5.1 Reverse flows 

A second solution of the F-S equation exists in the range J 3 m - n < /? < representing 
reverse flow. Figure 5 shows the variation of the velocity profile now for a 
negative shooting angle for /? = -0.01 and -0.1. As above, the profile passes 
through the asymptotic condition of unity, from below with increasing a. This 
behavior is again captured by a bisection procedure once the asymptotic limit is 
first surpassed. Table 7 gives a- values for reversed flow for the Rvalues 
considered by Katagiri in Table 6. This table required 1 10s of computing time and 
should be correct to several units in the last place. 

In the interest of adding an ending to a previous investigation on reverse flow, we 
consider one of the earliest and most influential articles regarding the shooting 
method [9]. In Table 8, the Cebeci and Keller investigation is completed. In Ref. 
9, only 3 digits for the given /Ts were found for reverse flow. Except for the last 3 
entries, the 3 digits are in nearly complete agreement with the present calculation. 
The last entry was modified since the value in Ref. 9 was in the region of no 
solution. It should be noted that the value close to the existence boundary on the 
lower branch nearly matches that of the upper branch as expected. Figure 6 
displays the velocity profiles for the cases of Ref. 9. The region of reversed flow is 
clearly evident. 

6. CONCLUSION 

An efficient algorithm, arguably the most efficient to date, has been devised to 
resolve the Falkner-Skan equation of boundary layer theory. The algorithm is 
simple in that it is based on a Maclaurin series solution through recurrence. Most 



Table 6 

Comparison to Katagiri Results 



JL 



1.00E+00 


1.23258765682E+00 


9.50E-01 


1.20546125458E+00 


9.00E-01 


1. 1777278 1917E+00 


8.50E-01 


1.14934554396E+00 


8.00E-01 


1.12026765738E+00 


7.50E-01 


1.090441 562 17E+00 


7.00E-01 


1.05980777320E+00 


6.50E-01 


1.02829859292E+00 


6.00E-01 


9.95836440616E-01 


5.50E-01 


9.62331717602E-01 


5.00E-01 


9.27680039837E-01 


4.50E-01 


8.91758591637E-01 


4.00E-01 


8.54421231 190E-01 


3.50E-01 


8.15491778664E-01 


3.00E-01 


7.7475458031 1E-01 


2.50E-01 


7.31940848513E-01 


2.00E-01 


6.86708181032E-01 


1.50E-01 


6.386085 12560E-01 


1.00E-01 


5.87035219198E-01 


5.00E-02 


5.31129630465E-01 


0.00E+00 


4.69599988361E-01 


-1.00E-02 


4.56454824718E-01 


-2.00E-02 


4.42982276168E-01 


-3.00E-02 


4.29156500469E-01 


-4.00E-02 


4.14947955734E-01 


-5.00E-02 


4.00322595418E-01 


-6.00E-02 


3.85240819783E-01 


-7.00E-02 


3.69656086478E-01 


-8.00E-02 


3.53513033028E-01 


-9.00E-02 


3.36744882009E-01 


-1.00E-01 


3.19269759843E-01 


-l.lUb-Ul 


3.0098531 1 136E-01 


-1.20E-01 


2.81760524240E-01 


-1.30E-01 


2.61422755987E-01 


-1.40E-01 


2.397359553 12E-01 


-1.50E-01 


2.16361405647E-01 


-1.60E-01 


1.90779855269E-01 


-1.70E-01 


1.621 14677053E-01 


-1.80E-01 


1.28636220596E-01 


-1.90E-01 


8.56997440597E-02 




Fig. 5. Variation of the velocity profile for reversed flow for 
/? = -0.01 and -0.1. 



Table 7 

Reversed flow for 0.198837735 < J3 < 






a 


-1.00E-02 


-4.23209260392E-02 


-2.00E-02 


-6.51685855429E-02 


-3.00E-02 


-8.25629651 188E-02 


-4.00E-02 


-9.66367874086E-02 


-5.00E-02 


-1.08271083286E-01 


-6.00E-02 


-1.17924111038E-01 


-7.00E-02 


-1.25858531722E-01 


-8.00E-02 


-1.32227654380E-01 


-9.00E-02 


-1.37113811802E-01 


-1.00E-01 


-1.40546212979E-01 


-1.10E-01 


-1.42507910309E-01 


1 A1 

-1.20E-01 


1 A O C\1 C 1 C\ A 1 C OT7 A1 

- 1 .42935 1 94358E-0 1 


-1.30E-01 


- 1.41 709502 152E-01 


-1.40E-01 


-1.38638925254E-01 


-1.50E-01 


-1.33421237895E-01 


-1.60E-01 


-1. 255676642 17E-01 


-1.70E-01 


-1.14227239098E-01 


-1.80E-01 


-9.7692060 1346E-02 


-1.90E-01 


-7.13359060034E-02 



importantly, the a highly accurate series evaluation is enabled through a Wynn- 
epsilon convergence acceleration coupled to continuous analytical continuation. 
The coupling has been demonstrated to provide the solution far into the asymptotic 
region. We then follow a shooting strategy based on bisection to force the 
shooting angle a (skin friction coefficient) toward its correct value as the velocity 
relaxes to its asymptotic profile. There is no Newton-Raphson search to adjust and 
control; and as long as the shooting angle is above its expected value, the initial 
guess is arbitrary. For an example of robustness, consider /3 = 1000. An initial 
value for the iteration of a = 200 gives a converged value of 36.5171968. The 
accuracy of this value however awaits the development of second algorithm of 
greater of equal accuracy—which will be a challenge if method developers 
continues to use differential equation solvers. The present algorithm also 
accommodates reversed flow an example of which is shown in Fig. 7 for /?= -0.16, 
where the F-S functions are contrasted for forward and reversed flow. 

The significance of the investigation presented here is that it closes the chapter on 
the development of highly accurate algorithms for the F-S equation. The next 
challenge is to develop the algorithm for the non-adiabatic case, including the 
energy equation. Also under consideration is the case for stretching surfaces with 



blowing and suction. Of particular note is that the Maclaurin series solution, as 
suggested by Blasius, has been vindicated as the solution representation of choice 
and may prove to be a viable approach for resolution of more comprehensive flow 
scenarios in the future. 



Table 8 

Completing calculations for Cebeci and Keller [9] 



investigation 








a 




-9.16200E- 


-03 


-3.99994660387E- 


-02 


-2.47890E- 


•02 


-7.39993700921E- 


-02 


-4.02860E- 


•02 


-9.70005720161E- 


-02 


-4.97450E- 


•02 


-1.08000456249E- 


-01 


-1.01763E- 


•01 


-1.41000029373E- 


-01 


-7.95960E- 


-02 


-1.31999450790E- 


-01 


-1.52118E- 


•01 


-1.31999333488E- 


-01 


-1.80553E- 


-01 


-9.65623555697E- 


-02 


-1.80552E- 


-01 


-9.65644238587E- 


-02 


-1.96348E- 


-01 


-4.00005870180E- 


-02 


-1.98826E- 


-01 


-2.88367895808E- 


-03 


-1.98837735E-01 


-5.77014 E 


-06 



1.2 



1.0 



0.8 - 



0.6 



J5- 
4- 



J3= -0.005913 



0.4 



0.2 - 



0.0 - 



-0.2 



J3= -0.198837735 




10 



P 

5.91300E-03 
7.46500E-03 
9.16200E-03 
1.30000E-02 
2.47890E-02 
4.02700E-02 
4.97450E-02 
7.95960E-02 
1.01763E-01 
1.52118E-01 
1.80552E-01 
1.98837735E- 

1 

15 



a 

-2.99983621 290E-02 
-3.5001 0528604E-02 
-3.99994660387E-02 
-4.99999989222E-02 
-7.39993700921 E-02 
-9.6980272221 8E-02 
-1.08000456249E-01 
-1.31999453439E-01 
-1.41000029373E-01 
-1.31999333488E-01 
-9.65644238587E-02 
01 -5.77014 E-03 

1 

20 



25 



7/ 



Fig. 6. Examples of reversed flow of Ref. 9. 
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